Code
library(tidyverse)
library(sf)
library(tmap)
library(terra)
library(stars)
library(kableExtra)December 3, 2025
Link to GitHub repository: https://github.com/aakriti-poudel-chhetri/impacts-of-extreme-weather
From February 13 to 17, 2021, Winter Storm Uri brought unusually cold temperatures, heavy snowfall, ice and freezing rain across Texas. Many places experienced their coldest conditions in more than three decades. The state’s power infrastructure was not designed for such extreme cold, therefore numerous power plants and grid components froze or malfunctioned. At the same time, electricity demand surged as residents attempted to heat their homes, placing severe pressure on the power system. To avoid a complete grid failure, authorities implemented rolling blackouts that were intended to last about an hour, but in reality many households lost electricity for several hours or even longer. An estimated 4.5 million homes experienced outages. These disruptions triggered additional problems, including widespread water and food shortages, burst pipes, and a lack of indoor heating. At least 246 deaths were linked directly or indirectly to the storm. Economic losses exceeded $195 billion, making Winter Storm Uri the expensive natural disaster in Texas history.
Understanding how communities were affected is crucial for preparing for future extreme weather events. It can help identify which neighborhoods are most vulnerable and guide decisions about where to focus recovery and emergency resources.
Thus, this study examines how Winter Storm Uri affected communities across the Houston metropolitan area by estimating the number of homes that lost power and assessing whether lower-income neighborhoods faced a disproportionate share of the impacts. Using satellite-based nighttime lights data from the VIIRS VNP46A1 product, we identify areas where electricity was likely lost by comparing light intensity before and after the storm. These blackout areas are then combined with OpenStreetMap building data to estimate the number of affected homes and linked with U.S. Census Bureau data to evaluate socioeconomic differences—such as variations in median household income—between impacted and non-impacted census tracts.
Understanding how power outages unfolded across Houston during Winter Storm Uri is important for evaluating both the extent of the blackout and the communities most heavily affected. Using satellite-based detection of outage areas alongside socioeconomic data, this study examines whether the storm’s impacts were distributed unevenly across neighborhoods. The main research question guiding this analysis is:
“How were power outages during Winter Storm Uri distributed across the Houston metropolitan area, and did these outages disproportionately affect lower-income neighborhoods?”
This study also addresses the following sub-questions:
For this analysis, data were extracted from the multiple sources.
VIIRS data is distributed through NASA’s Level-1 and Atmospheric Archive & Distribution System Distributed Active Archive Center (LAADS DAAC). Many NASA Earth data products are distributed in 10x10 degree tiles in sinusoidal equal-area projection. Tiles are identified by their horizontal and vertical position in the grid. Houston lies on the border of tiles h08v05 and h08v06. These two tiles per date were pre-downloaded and provided by the team.
To assess the impacts of Winter Storm Uri on Houston communities, a multi-step approach combining spatial and socioeconomic data was adopted. This allowed us to identify areas affected by power outages, estimate the number of homes impacted, and examine patterns across different socioeconomic groups. The major steps of the analysis are summarized below:
Let’s start by importing all the necessary libraries.
First, we read our datasets.
# Read night lights data
tile05_07 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021038.h08v05.001.2021039064328.tif'))
tile06_07 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021038.h08v06.001.2021039064329.tif'))
tile05_16 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather','data', 'VNP46A1', 'VNP46A1.A2021047.h08v05.001.2021048091106.tif'))
tile06_16 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021047.h08v06.001.2021048091105.tif'))
# Read roads data
highways <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'gis_osm_roads_free_1.gpkg'),
query = "SELECT * FROM gis_osm_roads_free_1 WHERE fclass = 'motorway'") %>%
st_transform(crs = 'EPSG:3083')
# Read houses data
houses <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data',
'gis_osm_buildings_a_free_1.gpkg'),
query = "SELECT *
FROM gis_osm_buildings_a_free_1
WHERE (type IS NULL AND name IS NULL)
OR type IN ('residential', 'apartments', 'house', 'static_caravan', 'detached')") %>%
st_transform('EPSG:3083')
# Explore the contents of the geodatabase
socioeconomic <- st_layers(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'))
# Extract the census tract information
census_tract <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'), layer = 'ACS_2019_5YR_TRACT_48_TEXAS') %>%
st_transform('EPSG:3083')
# Extract the income information
income <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'), layer = 'X19_INCOME')We will merge the raster tiles and crop it to the Houston bounding box to focus the analysis on the study area. It will ensure only valid night-light data within Houston were used for assessing power outages.
# Merge two raster data of February 7
nightlight_07 <- st_mosaic(tile05_07, tile06_07)
# Merge two raster data of February 16
nightlight_16 <- st_mosaic(tile05_16, tile06_16)
# Create Houston bounding box
houston_bbox <- st_bbox(c(xmin = -96.5, xmax = -94.5, ymin = 29, ymax = 30.5),
crs = 'EPSG:4326')
# Crop all raster data to Houston area
nightlight_07_crop <- st_crop(nightlight_07, houston_bbox)
nightlight_07_crop[(nightlight_07_crop > 1000) | (nightlight_07_crop <= 0)] <- NA
nightlight_16_crop <- st_crop(nightlight_16, houston_bbox)
nightlight_16_crop[(nightlight_16_crop > 1000) | (nightlight_16_crop <= 0)] <- NAHere, we visualize the comparison between night-light intensities in Houston before and after the February 2021 winter storm. By mapping the February 7 (pre-storm) and February 16 (post-storm) night-light data side by side, the analysis shows areas where lighting and likely power availability changed, helping to identify regions affected by the blackout.
tmap_options(component.autoscale = FALSE)
# Map for February 7 (before storm)
beforestorm <- tm_shape(nightlight_07_crop) +
tm_raster(col.scale = tm_scale_continuous(values = 'inferno'),
col.legend = tm_legend(show = TRUE,
title = 'Night light intensity before storm',
title.size = 1,
orientation = 'landscape',
width = 35,
height = 3,
legend.just = c(0, 1)))+
tm_title('Houston Night Lights\nFebruary 07, 2021 (Before storm)',
size = 1.2, fontface = 'bold') +
tm_compass(position = c('left', 'bottom'),
type = 'arrow',
size = 1.5,
text.size = 0.7,
color.dark = 'grey50',
text.color = 'white') +
tm_scalebar(position = c('left', 'bottom'),
text.size = 0.8,
color.dark = 'grey50',
text.color = 'white') +
tm_layout(bg.color = 'white',
outer.bg.color = 'white',
frame = TRUE,
legend.show = TRUE) +
tm_graticules(x.show = TRUE, y.show = FALSE, lwd = 0.2, col = "white", alpha = 0.3, labels.size = 0.5)
#tm_graticules(lwd = 0.2, col = "white", alpha = 0.3, labels.size = 0.5)
# Map for February 16 (after storm)
afterstorm <- tm_shape(nightlight_16_crop) +
tm_raster(col.scale = tm_scale_continuous(values = 'inferno'),
col.legend = tm_legend(show = TRUE,
title = 'Night light intensity after storm',
title.size = 1,
orientation = 'landscape',
width = 35,
height = 3)) +
tm_title('Houston Night Lights\nFebruary 16, 2021 (After storm)',
size = 1.2, fontface = 'bold') +
tm_compass(position = c('left', 'bottom'),
type = 'arrow',
size = 2.5,
text.size = 0.7,
color.dark = 'grey50',
text.color = 'white') +
tm_scalebar(position = c('left', 'bottom'),
text.size = 0.8,
color.dark = 'grey50',
text.color = 'white') +
tm_layout(bg.color = 'white',
outer.bg.color = 'white',
frame = TRUE,
legend.show = TRUE) +
tm_graticules(x.show = TRUE, y.show = FALSE, lwd = 0.2, col = "white", alpha = 0.3, labels.size = 0.5)
tmap_arrange(beforestorm, afterstorm, ncol = 2)We will now create a blackout mask by identifying areas of Houston where night light intensity decreased significantly after the storm, indicating potential power outages. By calculating the difference between before and after night light, cropping to the study area, filtering out minor changes, and vectorizing the result, the mask provides a spatial representation of blackout affected regions for further analysis.
# Calculate the difference in light intensity
nightlight_diff <- nightlight_07 - nightlight_16
# Crop and reclassify the raster data
mask <- st_crop(nightlight_diff, houston_bbox)
mask[mask < 200] <- NA
# Vectorize the blackout mask
blackout <- st_as_sf(mask,
as_points = FALSE,
merge = TRUE) %>%
st_make_valid() %>%
st_transform(crs = 'EPSG:3083')In this step, highway geometries are combined and buffered to exclude areas where changes in night light intensity could be due to traffic rather than power outages. By removing blackout areas within 200 meters of highways, the analysis isolates regions where lighting loss more likely reflects true residential or commercial power outages.
# Combine all highway geometries into one
highways_union <- st_union(highways)
# Create 200m buffer around all highways
highways_buffer <- st_buffer(highways_union, dist = 200)
# Find areas that experienced blackouts further than 200m from highways
blackout_far <- st_difference(blackout, highways_buffer)We visualize the spatial extent of power outages across Houston by mapping the blackout areas. The map highlights the neighborhoods affected, providing a clear view of the geographic distribution of the electricity disruption during the storm.
# Visualize the blackout in Houston
tm_basemap('OpenStreetMap') +
tm_shape(blackout_far) +
tm_polygons(fill = 'ivory', fill_alpha = 0.8, col = 'black') +
tm_title('Blackout areas in Houston', fontface = 'bold', size = 1) +
tm_layout(legend.show = TRUE) +
tm_scalebar(position = c('left', 'bottom')) +
tm_compass(position = c('right', 'top')) +
tm_graticules(lwd = 0.2, alpha = 0.3, labels.size = 0.5) +
tm_add_legend(labels = c("Blackout Areas"),
fill = c('black'),
type = 'polygons',
position = c('left', 'top'))We estimate the number of homes in Houston that lost power during the February 2021 winter storms. By identifying which building footprints intersect with the blackout mask, the analysis isolates the homes likely affected by outages. Converting buildings to centroids allows for precise spatial calculations and comparisons between impacted and non-impacted houses, resulting in an estimated total of approximately 157,970 homes without power.
# Keep buildings that intersect with blackout areas
houses_blackout <- st_intersects(houses, blackout_far)
# Check if each building is in a blackout area
in_blackout <- lengths(houses_blackout) > 0
# Filter buildings in blackout areas
blackout_houses <- houses[in_blackout, ]
# Number of impacted buildings
n_blackout_houses <- nrow(blackout_houses)
# Convert all buildings to centroids
houses_points <- houses %>%
st_centroid()
# Houses impacted by blackouts
houses_impacted <- blackout_houses %>%
st_centroid()
# Houses not impacted by blackouts
houses_not_impacted <- houses_points %>%
filter(!osm_id %in% blackout_houses$osm_id)Let’s create a table to show the number of houses affected by the power outage, making it easier to visualize the impact.
# Table of homes that experienced blackouts
houses_df <- houses_impacted |>
group_by(type) |>
summarise(count = n(),
percent = round((n() / nrow(houses_impacted)) * 100, digits = 2)) |>
st_drop_geometry()
# Add total row
houses_df <- houses_df |>
add_row(type = 'Total', count = sum(houses_df$count))
# Change type NA to character string
houses_df$type[6] <- "NA"
# Display kable table
options(knitr.kable.NA = "")
kbl(houses_df,
table.attr = 'class="custom-table"',
col.names = c(" Types", "Count", "Percentage (%)"),
align = 'c')| Types | Count | Percentage (%) |
|---|---|---|
| apartments | 1136 | 0.72 |
| detached | 353 | 0.22 |
| house | 19760 | 12.51 |
| residential | 1395 | 0.88 |
| static_caravan | 80 | 0.05 |
| NA | 135243 | 85.61 |
| Total | 157967 |
Plot map for the houses that were impacted by the power outage.
# Create the map
tm_basemap('OpenStreetMap') +
tm_shape(houses_not_impacted) +
tm_dots(fill = 'midnightblue', size = 0.02, fill_alpha = 0.8) +
tm_shape(houses_impacted) +
tm_dots(fill = 'firebrick', size = 0.02, fill_alpha = 0.8) +
tm_title('Houses in Houston that lost power',
position = tm_pos_out('center', 'top'),
fontface = 'bold', size = 1.2) +
tm_layout(legend.outside = TRUE,
legend.outside.position = 'right',
legend.title.size = 1,
legend.text.size = 0.8) +
tm_add_legend(labels = c('Impacted houses', 'Not impacted houses'),
fill = c('firebrick', 'midnightblue'),
col = c('firebrick', 'midnightblue'),
alpha = c(0.2, 0.7, 0.5),
title = 'Legend') +
tm_scalebar(position = c('left', 'bottom'),
text.size = 0.6) +
tm_compass(position = c('right', 'bottom'),
type = '4star',
size = 2)The following process links census tract boundaries with median household income data to analyze the socioeconomic characteristics of neighborhoods affected by the blackout. By joining the income dataset to the census tract geometries, the analysis enables comparisons between impacted and non-impacted areas and supports investigation of whether power outages disproportionately affected higher- or lower-income communities.
Let’s verify CRS. Verifying CRS alignment is essential for accurate spatial analysis, as mismatched projections can lead to incorrect geometry intersections, area calculations, or mapping errors.
Warning: Update coordinate reference systems to match!
By aligning the coordinate systems and cropping to the Houston study area, census tracts and blackout affected houses are prepared for spatial comparison. Intersections between the two datasets identify which tracts experienced power outages, and a new column flags these impacted areas. Filtering to only the affected tracts enables a targeted analysis and visualization of how the blackout influenced neighborhoods across the city.
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
# Transform and clean blackout houses
blackout_houses_transformed <- blackout_houses %>%
st_transform(crs = 'epsg:4326')
# Identify the blackout areas
census_blackout <- st_intersects(census_income_crop, blackout_houses_transformed)
# Create the column
census_blackout_col <- census_income_crop %>%
mutate(blackout = lengths(census_blackout) > 0)
# Select only the values that are true for plotting
census_blackout_col_filtered <- census_blackout_col[census_blackout_col$blackout, ]The following map highlights Houston census tracts that experienced power outages during the February 2021 winter storms. Affected tracts are shown in red, providing a clear visualization of the spatial distribution of blackouts across the city. This allows for easy identification of neighborhoods impacted by the storm and supports further analysis of socioeconomic patterns in relation to power loss.
# Map the census tracts that lost power
tm_basemap('CartoDB.VoyagerNoLabels') +
tm_shape(census_blackout_col_filtered) +
tm_polygons(fill = 'blackout',
fill.scale = tm_scale_categorical(
values = c('TRUE' = 'firebrick'),
labels = c('TRUE' = 'Blackout')),
fill.legend = tm_legend(title = 'Blackout status'),
col = 'grey50',
lwd = 0.2) +
tm_title('Census tracts in Houston that lost power',
position = tm_pos_out('center', 'top'),
fontface = 'bold',
size = 1) +
tm_layout(legend.outside = TRUE,
legend.outside.position = 'right',
legend.bg.color = 'white',
legend.bg.alpha = 0.8) +
tm_compass(position = c('right', 'bottom'),
type = 'arrow',
size = 1) +
tm_scalebar(position = c('left', 'top'),
text.size = 0.8)This plot compares the distributions of median household income between Houston census tracts that experienced blackouts and those that did not during the February 2021 winter storms. By visualizing income differences, it allows for assessment of whether power outages disproportionately affected higher- or lower-income neighborhoods and provides insight into potential socioeconomic disparities in storm impacts.
ggplot(census_blackout_col %>% filter(!is.na(B19013e1)), # Remove NA values
aes( x = blackout, y = B19013e1, fill = blackout)) +
geom_boxplot(alpha = 0.8) +
scale_x_discrete(labels = c('TRUE' = 'Blackout',
'FALSE' = 'No blackout')) +
scale_fill_manual(values = c('TRUE' = 'firebrick',
'FALSE' = 'midnightblue'),
labels = c('Blackout', 'No blackout')) +
labs(title = 'Median household income for census tracts',
subtitle = 'Blackout status in Houston census tracts',
x = 'Blackout status',
y = 'Median household income in $',
fill = 'Status') +
scale_y_continuous(labels = scales::dollar_format()) +
theme_minimal() +
theme(plot.title = element_text(face = 'bold', size = 16),
plot.subtitle = element_text(size = 14),
legend.position = 'bottom')This analysis used satellite data to compare median household incomes between areas affected and unaffected by the blackout. Interestingly, higher-income census tracts were more affected, as their median income was slightly higher than that of non-blackout areas. This result may reflect several limitations such as wealthier neighborhoods often have more outdoor lighting, census tract–level data may overlook neighborhood differences, and blackout areas could be misclassified due to highway lighting or cloud cover. The analysis may also have assumed all structures were residential, while vegetation near power lines and the 200-meter buffer used could further limit accuracy of true lighting.
Melaku, Fares, and Awal (2023); Zhou et al. (2024); (choi2021?)
@online{poudel2025,
author = {Poudel, Aakriti},
title = {Identifying the Impacts of Extreme Weather},
date = {2025-12-03},
url = {https://aakriti-poudel-chhetri.github.io/posts/2025-12-impacts-of-extreme-weather/},
langid = {en}
}